install.packages("survival", dependencies=TRUE)
install.packages("mstate", dependencies=TRUE)

rm(list=ls())
library(foreign)
library(survival)
library(mstate)


###########################################
######## Models In the Appendix############
###########################################

###### Model with South Africa (Table 5): ######
rm(list=ls())
#Make Sure Replication Data is in the Same Folder as This File
dat.sa <- read.dta("use_rep.dta")

mod.sa <- coxph(Surv(t, status) ~ 
		nuk7set1_1 + npt_eff_1+   lmid_past5yrs_1 +  neopatr_1+
		nuk7set1_2 + npt_eff_2+   lmid_past5yrs_2+ neopatr_2+
		nuk7set1_3 + npt_eff_3 +  lmid_past5yrs_3+ neopatr_3
				+ strata(trans), data = dat.sa, method = "breslow")
summary(mod.sa)


###### Replication of Table 6 in the Appendix ######
rm(list=ls())
#Make Sure Replication Data is in the Same Folder as This File
dat.way <- read.dta("useWAY1_rep.dta")#This is Jo+Gartke+ Singh and Way with cases of pursuit

dat.way$mid.lt.1 <- dat.way$lmid_past5yrs_1 * log(dat.way$t + 1)

mod.way <- coxph(Surv(t, status) ~ 
		nuk7set1_1 + npt_eff_1+   lmid_past5yrs_1 +  neopatr_1 +mid.lt.1+
		nuk7set1_2 + npt_eff_2+   lmid_past5yrs_2+ neopatr_2+
		nuk7set1_3 + npt_eff_3 +  lmid_past5yrs_3+ neopatr_3
				+ strata(trans), data = dat.way, method = "breslow")
summary(mod.way)


###### Replication of Table 7 in the Appendix ######
rm(list=ls())
#Make Sure Replication Data is in the Same Folder as This File
dat.wayex <- read.dta("useWAY1Ex_rep.dta")#This is Jo+Gartke+ Singh and Way with cases of exploration and pursuit

dat.wayex$mid.lt.1 <- dat.wayex$lmid_past5yrs_1 * log(dat.wayex$t + 1)

mod.wayex <- coxph(Surv(t, status) ~ 
		nuk7set1_1 + npt_eff_1+   lmid_past5yrs_1 +  neopatr_1 +mid.lt.1+
		nuk7set1_2 + npt_eff_2+   lmid_past5yrs_2+ neopatr_2+
		nuk7set1_3 + npt_eff_3 +  lmid_past5yrs_3+ neopatr_3
				+ strata(trans), data = dat.wayex, method = "breslow")
summary(mod.wayex)



#The following dataset is used to replicate Tables 8-16 of the Appendix
#Make Sure Replication Data is in the Same Folder as This File
rm(list=ls())
dat.1 <- read.dta("use_rep.dta")
dat.1 <- subset(dat.1, dat.1$ccode != 560)

#Control for security environment (Table 8):

mod.3 <- coxph(Surv(t, status) ~ 
		nuk7set1_1 + npt_eff_1+   lmid_past5yrs_1 +  neopatr_1 +nucriv_1+nucsecguart_1+
		nuk7set1_2 + npt_eff_2+   lmid_past5yrs_2+ neopatr_2+nucriv_2+nucsecguart_2+
		nuk7set1_3 + npt_eff_3 +  lmid_past5yrs_3+ neopatr_3+nucriv_3+nucsecguart_3
				+ strata(trans), data = dat.1, method = "breslow")
summary(mod.3)

# Control for Civil Nuclear Assistance (Table 9)
mod.4 <- coxph(Surv(t, status) ~ 
		nuk7set1_1 + npt_eff_1+   lmid_past5yrs_1 +  neopatr_1 + civass_1 + 
		nuk7set1_2 + npt_eff_2+   lmid_past5yrs_2+ neopatr_2+ civass_2 +
		nuk7set1_3 + npt_eff_3 +  lmid_past5yrs_3+ neopatr_3+ civass_3
				+ strata(trans), data = dat.1, method = "breslow")
summary(mod.4)

# Control for ENR (latency) (Table 10)
mod.6 <- coxph(Surv(t, status) ~ 
		nuk7set1_1 + npt_eff_1+   lmid_past5yrs_1 +  neopatr_1 + latency1_1+
	    nuk7set1_2 + npt_eff_2+   lmid_past5yrs_2+ neopatr_2+ latency1_2 +
		nuk7set1_3 + npt_eff_3 +  lmid_past5yrs_3+ neopatr_3 +latency1_3 
				+ strata(trans), data = dat.1, method = "breslow")
summary(mod.6)

#Control for Instability-1 Year Change (Table 11)
#Generate time interactions to correct for violations of the proportional hazards assumption
dat.1$mid.lt.1 <- dat.1$lmid_past5yrs_1 * log(dat.1$t + 1)
dat.1$polity1.1 <- dat.1$polity1yd_1* log(dat.1$t + 1)

mod.7 <- coxph(Surv(t, status) ~ 
		nuk7set1_1 + npt_eff_1+   lmid_past5yrs_1 +  neopatr_1 + polity1yd_1+ dat.1$mid.lt.1+dat.1$polity1.1+
		nuk7set1_2 + npt_eff_2+   lmid_past5yrs_2+ neopatr_2+ polity1yd_2+
		nuk7set1_3 + npt_eff_3 +  lmid_past5yrs_3+ neopatr_3+ polity1yd_3
				+ strata(trans), data = dat.1, method = "breslow")
summary(mod.7)

#Control for Instability-2 Year Change (Table 12)
#Generate time interactions to correct for violations of the proportional hazards assumption
dat.1$mid.lt.1 <- dat.1$lmid_past5yrs_1 * log(dat.1$t + 1)
dat.1$polity2.1 <- dat.1$polity2yd_1* log(dat.1$t + 1)

mod.8 <- coxph(Surv(t, status) ~ 
		nuk7set1_1 + npt_eff_1+   lmid_past5yrs_1 +  neopatr_1 + polity2yd_1+ dat.1$mid.lt.1+ dat.1$polity2.1+
		nuk7set1_2 + npt_eff_2+   lmid_past5yrs_2+ neopatr_2+ polity2yd_2+
		nuk7set1_3 + npt_eff_3 +  lmid_past5yrs_3+ neopatr_3+ polity2yd_3
				+ strata(trans), data = dat.1, method = "breslow")
summary(mod.8)


#Control for Instability-5 Year Change (Table 13)
#Generate time interactions to correct for violations of the proportional hazards assumption
dat.1$mid.lt.1 <- dat.1$lmid_past5yrs_1 * log(dat.1$t + 1)
dat.1$polity5.1 <- dat.1$polity5yd_1* log(dat.1$t + 1)

mod.9 <- coxph(Surv(t, status) ~ 
		nuk7set1_1 + npt_eff_1+   lmid_past5yrs_1 +  neopatr_1 + polity5yd_1+ polity5.1+
		nuk7set1_2 + npt_eff_2+   lmid_past5yrs_2+ neopatr_2+ polity5yd_2+
		nuk7set1_3 + npt_eff_3 +  lmid_past5yrs_3+ neopatr_3+ polity5yd_3
				+ strata(trans), data = dat.1, method = "breslow")
summary(mod.9)

#Control for global tech environment (Table 14)

mod.10 <- coxph(Surv(t, status) ~ 
		nuk7set1_1 + npt_eff_1+   lmid_past5yrs_1 +  neopatr_1 + nuk7set1_glsd_1 +
		nuk7set1_2 + npt_eff_2+   lmid_past5yrs_2+ neopatr_2+ nuk7set1_glsd_2+
		nuk7set1_3 + npt_eff_3 +  lmid_past5yrs_3+ neopatr_3+ nuk7set1_glsd_3
				+ strata(trans), data = dat.1, method = "breslow")
summary(mod.10)

# Control for Pressure (Table 15)
mod.11 <- coxph(Surv(t, status) ~ 
		nuk7set1_1 + npt_eff_1+   lmid_past5yrs_1 +  neopatr_1+
		nuk7set1_2 + npt_eff_2+   lmid_past5yrs_2+ neopatr_2+ pressure_2+
		nuk7set1_3 + npt_eff_3+  lmid_past5yrs_3+ neopatr_3+pressure_3
				+ strata(trans), data = dat.1, method = "breslow")
summary(mod.11)

#Control for NPT dummy instead (Table 16)
mod.12 <- coxph(Surv(t, status) ~ 
		nuk7set1_1 + npt_dum_1+   lmid_past5yrs_1 +  neopatr_1+
		nuk7set1_2 + npt_dum_2+   lmid_past5yrs_2+ neopatr_2+
		nuk7set1_3 + npt_dum_3 +  lmid_past5yrs_3+ neopatr_3
				+ strata(trans), data = dat.1, method = "breslow")
summary(mod.12)






#########################################
######## Time Varying Data Test #########
#########################################

rm(list=ls())
#Make Sure Replication Data is in the Same Folder as This File

dat.tv <- read.dta("usetv_npt_REP.dta")

#Generate time interactions to correct for violations of the proportional hazards assumption
dat.tv$mid.lt.1 <- dat.tv$lmid_past5yrs_1 * log(dat.tv$t + 1)
dat.tv$npt.lt.1 <- dat.tv$npt_effC_1 * log(dat.tv$t + 1)

dat.tv$mid.lt.2 <- dat.tv$lmid_past5yrs_2 * log(dat.tv$t + 1)
dat.tv$nuk.lt.2 <- dat.tv$nuk7set1_2 * log(dat.tv$t + 1)
dat.tv$npt.lt.2 <- dat.tv$npt_effC_2* log(dat.tv$t + 1)
dat.tv$neo.lt.2 <- dat.tv$neopatr_2 * log(dat.tv$t + 1)

#Time-Varying Multi-State Model (Table 17):
mod.tv <- coxph(Surv(t0, t, status) ~ 
nuk7set1_1 + npt_effC_1 + lmid_past5yrs_1 + neopatr_1	+ mid.lt.1+
nuk7set1_2 + npt_effC_2 + lmid_past5yrs_2 + neopatr_2 	+ mid.lt.2 + npt.lt.2 + nuk.lt.2+
nuk7set1_3  +  npt_effC_3 + lmid_past5yrs_3 + neopatr_3 +
strata(trans) + cluster(ccode), data = dat.tv, method = "efron")	
summary(mod.tv)
